Exact and Efficient Algoritiini to Discover Extreme Stochastic Events 
in Wind Generation over Transmission Power Grids 

Michael Chertkov, Mikhail Stepanov, Feng Pan, and Ross Baldick 



Abstract — In this manuscript we continue the thread of 
[M. Chertkov, F. Pan, M. Stepanov, Predicting Failures in 
Power Grids: The Case of Static Overloads, IEEE Smart Grid 
2011] and suggest a new algorithm discovering most probable 
extreme stochastic events in static power grids associated 
with intermittent generation of wind turbines. The algorithm 
becomes EXACT and EFFICIENT (polynomial) in the case of 
the proportional (or other low parametric) control of standard 
generation, and log-concave probability distribution of the 
renewable generation, assumed known from the wind forecast. 
We illustrate the algorithm's ability to discover problematic 
extreme events on the example of the IEEE RTS-96 model of 
transmission with additions of 10%, 20% and 30% of renewable 
generation. We observe that the probability of failure may grow 
but it may also decrease with increase in renewable penetration, 
if the latter is sufficiently diversified and distributed. 

I. Introduction 

In progress to becoming smarter, the power grid of today 
is undergoing multiple transformations [1], [2]. One of the 
envisioned changes consists in replacing a significant portion 
of the traditional fossil thermal plants by renewable gener- 
ation [3], in particular by solar and wind farms. This task, 
motivated by ecological and political reasons, will not be a 
simple substitution. The renewable sources of energy, will 
also be much less predictable and thus much more difficult 
to control. Both wind and solar fluctuate temporally and 
spatially. Even when forecasted, the renewable generation 
can be described only in probabilistic terms, suggesting that 
the existing toolbox of power engineering, which is largely 
deterministic, needs to be upgraded with computationally 
more challenging probabilistic tools. For any configuration 
of aggregated loads, which stays roughly constant for tens 
of seconds or even minutes, one ought to consider an 
ensemble of possible configurations of renewable generation 
undergoing significant and unpredictable changes during the 
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same time. This uncertain (but probabilistically predictable) 
configuration of renewable generation will be complemented 
by standard adjustable generation. The adjustment is required 
to stabilize the system and complete matching between re- 
quested load and total generation. Many configurations from 
this ensemble of matchings will be feasible, i.e. not violating 
transmission constraints. However, there will always be some 
number of probabilistically rare but strong fluctuations of the 
renewable sources for which an instantaneous adjustment 
of generation to loads will be problematic, as resulting in 
overload of some number of lines. Discovering extreme 
stochastic events of this rare but damaging type becomes an 
important practical problem which requires fast and largely 
non-existent algorithmic solutions. 

This manuscript is devoted to resolving (at least some part 
of) the aforementioned challenge. We assume that configu- 
ration of loads over the network is fixed, while generation is 
split into two parts: (a) the bigger, consisting of conventional 
generation, controlled and adjustable as a group, for example 
via common re-scaling (roughly representing, for example, 
so-called "droop" speed control with the proportion of each 
generator response to the frequency variation fixed, or rep- 
resenting "regulation" response when deployed regulation is 
also in proportion to the frequency variation); and (b) renew- 
able generation, say coming from wind farms, uncertain and 
not controlled, described in terms of a probability distribution 
function, assumed known from a forecast. We are interested 
to discover the most probable configuration(s) of the wind 
(more generally any fluctuating renewable source) which 
is troublesome, in terms of possibly violating transmission 
and/or controlled generation constraints, when no additional 
control efforts, such as curtailment of wind generation, load 
shedding, line switching, etc, are in place. With diversity 
of wind production across multiple wind sites, the extreme 
stochastic events will not be frequent, thus making the 
problem of finding the rare configuration causing the trouble 
as challenging as finding a needle in a haystack. To discover 
the extreme stochastic events we exploit and develop further 
the approach, originally introduced in theoretical physics, see 
e.g. [4], [5], [6], [7] then used to analyze performance of 
error-correction codes [8], and recently applied to predict 
static failures in power grids associated with fluctuations 
in loads [9]. We call this rare but most probable troubled 
instance of the wind the instanton. In a significant technical 
improvement in comparison with [9] this manuscript sug- 
gests a direct way, which is Exact and Efficient (polynomial 
in the size of the grid), to explore the structure of the problem 
for finding the instanton(s). This algorithmic improvement is 



achieved via mapping the instanton problem to minimization 
of a convex function, characterizing forecasted distribution of 
the wind, over the exterior of a tractable polytope, describing 
feasible solutions of the power flow equations. 

The material in the manuscript is organized as follows. 
We give some technical background, formulate the problem 
and briefly describe history of the instanton methodology 
in Section [III Section [1111 formulates the main theoretical 
result of the manuscript: exact and efficient algorithm for 
finding most probable wind extreme stochastic event for 
given configuration of the grid. The performance of the 



algorithm is illustrated in Section |IV] on example of the 
IEEE RTS96 system with 10%, 20% and 30% of renewable 
penetration. The results are summarized and discussed in 
Section |Vj where we also discuss future research challenges. 

II. Quasi-Static Power Flows, Wind Modeling 
AND Formulating the Instanton Problem 

Modern power grid is always in motion. Moreover many 
of the changes which take place on multiple time scales 
are inherently uncertain, being probabilistic in nature. With 
sufficient penetration of the renewable sources, wind turbines 
bring a particularly important source of fluctuations into the 
modern power grid. Production of wind can be forecasted 
only to a degree. Even with a perfect forecast, and due 
to a turbulent nature of the wind, one ought to describe 
output of the wind turbine integrated directly to the power 
grid (without curtailing) in probabilistic, rather than deter- 
ministic, terms. Relevant times scales, where fluctuations of 
wind generation dominate other sources of uncertainty (in 
particular these associated with fluctuations of demands), 
are in the range from tens of seconds to tens of minutes. 
At these temporal scales transient (sub-second) phenomena 
are already settled and quasi- static description of the power 
flows is appropriate [10], [11]. Taking a standard Direct 
Current (DC) power flow approximation (well justified for 
transmission network, where resistivity of power lines is 
significantly smaller that respective inductance, voltage vari- 
ations are small and change in phase between neighboring 
busses is small too) one arrives at the following relations 
between the vector of injected/consumed real powers, P, and 
the vector of phases, (p, over the power grid defined as a 
graph, ^ = (^0, ^i) (^0 and represent the set of nodes 
and the set of edges respectively), written in the form of 
conditions 



COND 



flow ' 



B(p=P, & (po = 



0, 

B={ -l/xij, {iJ}^Q\ 

^{hk}^Ql -1 

( -di, Qd \ 

Pn i ^ Qr 

V 0, ^jU^.U^g J 



(1) 

(2) 
(3) 



In the remainder of the paragraph we detail notations, and 
underlying notions, assumed in Eq. ([3]). 5 in Eq. ([T]) is 
the matrix of the graph Laplacian constructed from line 
inductive reactances, x, and the second condition in Eq. ([T]) 
sets the phase at an arbitrarily chosen (zero) node to zero. 
Different components of the vector, P, associated with nodes 
of the graph, are split in Eq. ^ into four groups: set of 
demand/consumption nodes, ^j; set of renewable nodes, 
set of standard generation nodes, and the remaining set 
of other (junction) nodes. Any component of the demand 
is assumed negative and fixed (on the time scale of interest 
measured in seconds-to-minutes). The renewable generation 
fluctuates according to the forecasted probability distribution 



lP(p)^exp(-5(p)), 



(4) 



where 5(p) is a known convex function of its multi- 
dimensional argument, achieving minimum at the most prob- 
able, equilibrium, configuration of the wind generation, Pq. 
This distribution function represents statistics of the wind 
fluctuations, collected over the time interval when loads 
do not change or change very little. (The assumption of 
5 convexity is reasonable for a sufficiently wide, and thus 
most probable, vicinity of the maximum output, in particular 
for popular modeling of the wind statistics via the Weibull 
function [12].) Standard (controllable) generation adjusts to 
match the imbalance between total consumption and total 
renewable generation (assuming that the loads stay constant). 
The adjustment is relatively fast (instantaneous within the 
quasi-static description) and we will simplify by assuming 
that the mismatch is shared between the controllable gener- 
ators in some pre-defined, automatic fashion. In particular, 
we will consider the case where the adjustment for each 
generator is proportional to its nominal generation, so that a 
entering Eq. ([3]), is 
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and it is the only degree of freedom on the standard gener- 
ation side which reacts absorbing changes in the renewable 
generation, p = (p/|/ G ^r)- We assume that a is set to 
unity at the equilibrium configuration, Pg. Then, {pi\i G ^g) 
constitutes output of the controllable generation preset in the 
beginning of the time interval of interest. Normally, this 
redispatch of the controllable generation is the output of 
the Optimum Power Flow (OPE) analysis, accounting for 
diversity in the generation cost at different sites, and executed 
periodically on the scale ranging from minutes to tens of 
minutes. The model of Eq. ^ schematically represents either 
droop control or regulation response, although details of both 
of these control modes differ from our model in details. 

One expects that a feasible Power Flow (PF) solution, 
representing a perturbed OPE, satisfies, in addition to the 
basic power flow relations ([T]), the following set of trans- 
mission (thermal) conditions representing line constraints on 
the amount of power which can flow safely trough the lines 



(without overheating or damage) 
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where Uij is the Hne {ij} rating. One also assumes that all 
the generators included in the proportional control are within 
their capacity bounds, which translates into the following 
cumulative constraints on the proportional control coefficient 
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where a and a are defined in accordance with the maxi- 
mum low and minimum high constraints, respectively, over 
individual controllable generators. In fact, typical control 
modes would allow for some generators reaching their limits; 
however, we simplify that issue here. 

Even though the most probable configuration of the 
wind generation corresponds to the power flow solution of 
COND flow, which is safely within the feasibility region of 
CONDedge^ COND power. Other less probable configurations 
of wind can violate one or more of the feasibility constraints. 
Naturally, we are first of all interested to discover the most 
probable instance of the wind, the instanton, which lies 
outside of the feasibility region: 



^inst 



: argminp5(p) ^^^^^ , (8) 
^int = Projection (COA^Z)((p,p,a))p , (9) 

COND = CONDfiowl^CONDedge^CONDpoweri^O) 



where Projection (COA^Z)((p,p, a) )p is the projection of the 
polytope COND to the p-space. Therefore, by construction 
(Dint is also a polytope. Alternatively we can rewrite Eq. ^ 
as 



argminp5(p) 



(11) 



where ^ext is a non-convex set, defined as an exterior of th^ 
convex set (polytope) "Dint, i.e. "Dejt = R^f''^ "Dint . Eq. (11 



states succinctly the instanton problem addressed in this 
manuscript. 

In more general formulations, not yet utilizing special 
structure of the optimization domain (D^xt specific to our 
problem, the instanton problem ( pTj ) can be solved within 
the machinery of non-convex optimization methods. For 
example, and as described in details in [13], [14], [9], one 
may search for the minimum of 5(p) with the help of 
a general-purpose optimization technique, specifically the 
downhill simplex (or "amoeba") method [15], [16]. Since 
the optimization domain is not concave, we expect to find 
many (candidate) instanton solutions. Each initialization of 
the instanton-amoeba could lead to a new instanton. The 
initialization selects a simplex, built on A/j + 1 feasible 
points (where A/j is the dimensionality of the p space, 
which is equal in our case to the number of renewable 
generators) from the error- surface separating (Dint and (Dext- 
Then the instanton-amoeba method evolves the simplex, via a 



sequential set of shifts, contractions and extensions, towards 
its eventual collapse to a local minimum of 5(p). Differ- 
ent random initiations will sample the space of instantons, 
thus generating the so-called instanton spectrum describing 
the frequency of a given instanton occurrence and also 
suggesting an estimation for the ordered list (with respect 
to their probability of occurrence and frequency) of top 
instantons. Repeated infinite number of times, the sampling 
would output the most probable instanton. However given 
that the number of initiations will be finite in reality, the 
most probable one (of the finite number of instantons found) 
gives a heuristic estimate from below for the probability of 
the most probable instanton. To ensure sampling quality one 
needs to continue random initiations till the most probable in- 
stantons would appear multiple number of times. (Typically, 
this require hundreds of initiations for the network measured 
in hundreds of nodes, of the type discussed below in Section 

lYl) 

However, the specific structure of our problem ( pT| ) allows 
a much faster and moreover exact resolution, than the one 
provided by the general purpose but computationally heavy 
instanton-amoeba method of [13], [14], [9]. As shown in the 
next Section, very specific features of ^ext, associated with 
the linear nature of the DC power flow and also with the 
single-parametric and linear control of fluctuations on the 
standard generation side, allow to solve Eq. ( pT) efficiently. 
The new approach, describe below, offers a significant 
algorithmic improvement in comparison with the general 
approach of the instanton-amoeba type. 

III. Exact and Efficient Algorithm to Discover 

THE INSTANTON(S) 

Given Eqs. ([T]) and Eq. ([5]), one can express (p and a via 
p, thus arriving at the following explicit (tractable polytope) 
expression for (Dint'. 
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where matrix B is the quasi-inverse of B accounting for 
the (po = constraint and thus regularizing the only zero 
eigenvalue of B. 

With this explicit formulation of (Dint the instanton prob- 
lem becomes tractable, as it reduces to a minimum over 
a tractable set of convex problems: 



min 



min5(p) 
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where K= \CONDedge^CONDpower\', Ca stands for any of 
the inequality constraints in COND ^^^^ U COND p^^^,^/, and 
is the saturated version of Ca by replacing inequality 



with equality. When the feasibihty set of a sub-problem 
a in Eq. (15) is empty we formally set to infinity. 
Eq. ( p?) is computationally tractable because it splits into 
K convex optimizations. The proof of the transition from 



Eq. ( pJ) to Eq. ( p3] ) is straightforward and it simply relays 
on testing the faces of (Dint sequentially. (See, e.g., detailed 
discussion of similar problem in [17].) For any of the internal 
minimizations in Eq. ( p3] ), which is feasible, having one facet 
of (Dint saturated is guaranteed (by construction) and unless 
there is a degeneracy in 5(p) there will be only one saturated 
facet per any minimization problem f^. 

Two comments are in order. First, our scheme is general, 
and can account for other types of the generation control 
besides adjusting generation linearly in response to the 
variation in p. However, one should also expect that the 
difficulty of the generalized problem will grow exponentially 
with the number of the control degrees of freedom. This 
exponential explosion is associated with the fact that pro- 
jection of the respective generalization of COA^Z)((p,p,a), 
where a is now multi-parametric, will generate generalized 
(Dint, which is much more complex as described via an 
exponentially large (in the dimensionality of a) number of 
contraints/inequalities. This is in spite of the fact that the 
original polytope, COA^Z)((p,p,a), is tractable. (Indeed, it is 
well known that projection of a tractable polytope along a 
subspace results in a polytope characterized, in the worst 
case, in terms of the set of constraints which is exponentially 
large in the size of the subspace. See [18], [19] for examples 
and related algorithmic discussions.) Second, given that the 
instanton problem is split into K tractable optimizations in 
Eq. ( p3) ), one also finds not only the instanton itself (absolute 
minimum) but also the ranked list of other extreme stochastic 
events associated with other faces of (Dint • We will call these 
other instantons, second-, third- etc according to their ranking 
in the derived hierarchy. This ranking is useful to discover the 
list of extreme stochastic events. Note, however, that strictly 
speaking this ranking does not necessarily corresponds to the 
actual ranked list of all possible extreme stochastic events. 
There are two reasons for that. First, any point from a small 
continuous vicinity of the (top ranked) instanton, sitting at 
the same face of (Dint as the instanton, will have a lower 5(p) 
weight than other instantons. (This is under assumption that 
the situation is not degenerate and the instanton corresponds 
to an interior point of the face.) Second, one may find a 
(discrete) configuration, which will be second in ranking 
within a minimization ranked k\ in Eq. (p?), but will still 



have a lower 5(p) weight than the top result of another 
minimization ranked fe, even when k\ < fe. 

IV. Numerical Example 

We test the algorithm on the standard IEEE RTS-96 model 
[20] extended with renewables. The results are shown in 
Figs. 1|2|3|4[ ) and commented upon below. To imitate effect 
of renewables, we took the base configuration of the model 
(standard generators, loads and transformer nodes are kept 
according to the data from [20]). Then, we add 3,6 and 9 
additional nodes, for new wind generators and connect each 
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Fig. L Graph of the IEEE RTS-96 model with nine new renewable 
generators (#73-#81) added. 
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(a) Configuration of renewable generation for the top 
three instantons. 
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(b) Top instantons are shown. Vertexes marked in 
red/orange/blue colors are of the most stressed renew- 
able generator from the instanton output. Edges marked 
red/orange/blue colors are the saturated ones for the 
respective instantons. 

Fig. 2. Instanton(s) for base configuration -i- 3 renewable generators: 10% 
of renewable penetration in average production. The instantons are ordered 
according to their cost values, *S(p). 
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(a) Configuration of renewable generation for the top 
three instantons. 




(b) Top instantons are shown. Vertexes marked in 
red/orange/blue colors are of the most stressed renew- 
able generator from the instanton output. (Bus #77 is the 
most stressed in the instanton #1 and the instanton #3.) 
Edges marked red/orange/blue colors are the saturated 
ones in the respective instantons. 

Fig. 3. Instanton(s) for base configuration + 6 renewable generators: 20% 
of renewable penetration in average production. The instantons are ordered 
according to their cost values, S{p). 
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(a) Configuration of renewable generation for the top 
three instantons. 




(b) Top instantons are shown. Vertexes marked in 
red/orange/blue colors are of the most stressed renew- 
able generator from the instanton output. (Bus #77 is 
the most stressed in the instanton #1 and instanton #3.) 
Edges marked red/orange/blue colors are the saturated 
ones in the respective instantons. 

Fig. 4. Instanton(s) for base configuration -i- 9 renewable generators: 30% 
of renewable penetration in average production. The instantons are ordered 
according to their cost values, S{p). 



of them to some number of other randomly selected nodes, 
where number of connections per new node is distributed 
according to the degree distribution of generators in the base 
configuration. To facilitate comparison, we also create the 
three configurations sequentially, such that the bigger one 
is built on the top of the smaller one. The new graph is 
shown in Fig. ([T]). The capacities of added lines are chosen 
equal (and roughly correspond to the median capacity of the 
existing lines). Inductances of the new lines are distributed 
uniformly in a median range of inductances of the base case. 
We choose a = and a = 2. We pick the simplest possible 
model for statistics of renewable generation, assuming that 
p is site-uncorrelated, positive (component by component) 
Gaussian, thus represented by 
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where one chooses p = (p^ | / G in a way that the total of 
typical renewable generation, L/e^i^Pr corresponds to 10%, 
20% and 30% (for the three test cases respectively) of the 
standard generation of the RTS-96 base case. Selecting p we 
also make sure that the resulting configurations are all in the 



regime where no transmission or generation constrains are 
violated, i.e. p G (Dint. Note, that the choice of the positive 
Gaussian and site-uncorrelated distribution is made here for 
illustrative purposes only. Actual correlations of wind will 
be more elaborate and interesting. However, one expects that 
the realistic, 5(p), representing actual wind forecast, will still 
be a convex function of p, thus making application of our 
algorithm to the more realistic situation as straightforward 
as for the synthetic positive Gaussian case discussed here. 

The results of our numerical tests, which main goal was to 
illustrate utility of the algorithm as of an exact and fast tool, 
are shown in Figs. ( 2|3|4| ) for the three levels of the renewable 
penetration respectively. For each example we present two 
figures showing in, (a) configuration of renewable generation 
p for the top three leading instantons in comparison with the 
base case, and in (b) structure of the top three instantons 
(renewable sites with the largest values of p/p marked) and 
respective saturated edges. A complete instanton analysis of 
an instance (graph+distribution) is very fast, it takes few 
seconds on a laptop. 

Here is the summary and brief discussion of the simulation 
results: 



• In the cases correspondent to 10% and 30% of renew- 
able penetration the top instantons are well localized 
on a site/node, in a sense that production at this single 
site is significantly larger than the typical value while 
deviation in production at the other renewable sites are 
much weaker. On the contrary, the top instanton is of a 
de-localized type in the intermediate case of 20% of 
renewable penetration. Note that the cost of the top 
instanton in the intermediate 20% case , S{p) = 1.04, is 
significantly lower than in the 10% case, S{p) = 11.28 
and also smaller than in the 30% case, S{p) = 2.49. 
This non-monotonicity of the cost with increase of 
the renewable penetration illustrates the generally very 
important point that addition of renewables can stress 
(transition from the 10% case to the (30% case) or de- 
stress (transition from the 10% case to 20% case) the 
network in terms of stochastic transmission congestion, 
depending on how it is done. 

• We also monitor the value of a at the instantons, and 
observe that it decreases with increase of the level of 
renewables: 0.75, 0.63 and 0.54 for the top instantons 
in the 10%, 20% and 30% cases respectively. This ex- 
presses the fact that the addition of renewables translates 
into reduction of the standard generation. 

• There is always one saturated edge per instanton. 
Somehow remarkably, this saturated edge is always 
one of the "old" edges (present in the base structure), 
and it is also positioned in majority of cases relatively 
far from the renewable sites which over-generate at 
the instanton. This observation emphasizes non-locality 
(and thus intrinsic difficulty) of the problem. 

• The standard generation constraints are not violated at 
the instantons, i.e. for each of the instantons found, 
a = 0<a<a = 2. a is unity at the base case and 
it (naturally) decreases with the level of renewable 
generation increase. 

V. Conclusions and Path Forward 
Summarizing, the main results of this work are: 

• We posed the problem of discovering extreme stochastic 
events in power grids, the instanton, associated with 
transmission overflows caused by fluctuations of renew- 
able generation. 

• We showed that the aforementioned problem is com- 
putationally tractable within the DC power flow setting 
and with a low parametric linear control of standard 
generation, for example of the proportional type. 

• We illustrated algorithmic utility and efficiency of our 
newly suggested algorithm on example of the IEEE 
RTS-96 grid with 10%, 20% and 30% of added re- 
newable generation. Main conclusions of our numerical 
tests are: (a) Addition of renewables may lead to signif- 
icant increase of the probability of failures (destructive 

^Note that this observation does not violate the N-l contingency 
requirement (a security measure enforcing that the system has a feasible 
solution for any one of the N edges of the graph removed) simply because 
the condition is enforced only at the equilibrium point. 



effect), but it may also help to make the network to 
become prone to failures (constructive effect), depend- 
ing on quantitative details of how the grid extension is 
done, (b) The instanton configurations represent global 
correlations in the graph, which shows itself in the fact 
that the (single) overloaded line is typically relatively 
far from the over-producing renewable sites, (c) Our 
algorithm is very efficient computationally in providing 
the exact assessment of the gain (or loss) associated 
with addition of renewables. 

We plan to continue this work on discovering extreme 
stochastic events in the power grids efficiently. Our main goal 
here is to design a reliable predictive tool, extreme stochastic 
event toolbox, capable to provide awareness (fast guidance 
to utility operator) in terms of predicting dangerous extreme 
stochastic events associated with the renewable generation. 
More specifically, this work will be continued along the 
following lines. 

• First of all, we will be testing scalability of the approach 
for larger, continental scale, transmission networks. Our 
task here is to approach complexity which would scale 
linearly with the size of the system. (This may be 
achieved, e.g. by replacing standard convex optimization 
solvers, by their linear scaling and distributed proxies.) 

• Our modeling of wind statistics needs to be more 
realistic. We plan to apply our algorithm to wind data 
taken from a realistic forecast, e.g. of the type available 
at [21], in particular accounting for long spatial (usu- 
ally, at least hundreds of kilometers long) correlations 
between different sites. Our future test beds will include 
models of ERCOT (with actual or planned wind farms 
in Western Texas), as well as models of other wind-rich 
parts of US power grids. 

• We will incorporate into the (so far static) scheme 
dynamic effects and actual (time series) measurements 
of the wind intensity. This will require extending the 
instanton approach to account for temporal Lagrangian 
correlations in the cost function, integrated in time 
over pre-history, and enforcing the transmission and 
generation conditions not only instantaneously but also 
over the (dicretized) time horizon. 

• Our model of proportional control is convenient for 
analysis, but it matches the details of droop control 
and of regulation response only schematically. We will 
modify the model to more fully represent these actions, 
in particular accounting for generation limits. 

• We envision incorporating the instanton analysis into 
control schemes, in particular in the tertiary (balance) 
control, for example penalizing top instanton configu- 
rations (and their vicinities) in the modified optimum 
power flows. Another interesting control option is to 
fit the instanton framework into an adversarial process, 
whereby any given control action can be matched by a 
corresponding instanton. This yields an iterative process 
(with each iteration roughly like the problem solved in 
the paper), which produces robust control actions and 



significant events, at the same time. 
• We will also work on extending the algorithm to the 
case of AC power flows (which is especially important 
in terms of accounting for additional issues related to 
voltage variations), e.g. utilizing new advances in related 
optimization techniques [22]. We will consider other 
schemes of transmission and voltage control, in particu- 
lar related to demand response [23], line switching [24], 
controlling DC tie-lines, capacitor banks, phase shifters 
and related [10]. We also plan to account for other type 
of stochastic issues, e.g. related to dynamic stability [25] 
and voltage collapse [26]. 
Finally, once the extreme stochastic event toolbox is devel- 
oped, we envision using it (as a black box) for developing 
new planning and control schemes for smart grids of the fu- 
ture, for example in the spirit of the general approach of [27] . 
We envision using this toolbox to solve related problems of 
discovering interdiction attacks on power grids [28], [29], 
and analyzing, controlling and preventing cascading failures 
[30], [31]. 
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